######################################################################
# Time-stamp: <05/07/06 14:59:04 ostolop>
######################################################################
package EP::EPC::Core::Data_Transformation::exp;

use strict;

use EP::Config;
use EP::Common;
use EP::Common::Error ':try';

use EP::Core::Data::Parser::exp::user;
use EP::Common::DBI::Dataset_Metadata_Type;

use PDL;
use PDL::NiceSlice;
use PDL::IO::FlexRaw;

use EP::Rclient;
use File::Temp;

# sub: process_transformation
#   Entry point for transformation processing.
# returns:
#   The new resulting dataset.
sub process_transformation {
  my $self = shift;
  my ( $dataset, $query ) = @_;

  if ( $query->param('transformation_intensity_to_ratio' ) ) {
    intensity_to_ratio ( $self,
			 $dataset,
			 $query->param('columns_per_exp'),
			 $query->param('channel1'),
			 $query->param('channel2'),
			 $query->param('log'),
			 $query->param('dye_swap') );
  } elsif ( $query->param('transformation_ratio_to_logratio' ) ) {
    ratio_to_logratio ( $self,
			$dataset,
			$query->param('log') );
  } elsif ( $query->param('transformation_avg_row' ) ) {
    avg_row_id ( $self, $dataset );
  } elsif ( $query->param ('transformation_impute_data_knn' ) ) {
    impute_data_knn ( $self, $dataset, $query->param('neighbours') );
  } elsif ( $query->param ('transformation_impute_data_replace_with_zeroes' ) ) {
    impute_data_replace_with_zeroes ( $self, $dataset );
  } elsif ( $query->param ('transformation_impute_data_replace_with_rowavgs' ) ) {
    impute_data_replace_with_rowavgs ( $self, $dataset );
  } elsif ( $query->param ('transformation_transpose' ) ) {
    transpose_data ( $self, $dataset );
  } elsif ( $query->param ( 'transformation_affy_rel_ratio' ) ) {
    affy_rel_ratio ( $self,
		     $dataset,
		     $query->param('log'),
		     $query->param('rel_to'),
		     $query->param('rel_to_column') );
  } elsif (  $query->param ( 'transformation_mean_center' ) ) {
    mean_center ( $self,
		  $dataset,
		  $query->param ( 'mean_center' ) );
  }

  delete $dataset->{data} and throw EP::Common::Error ( -value => CUD_EMPTY_SELECTION )
    if $dataset->data->matrix->isempty;

  # Copy the dataset now. Note: this causes $trf_dataset to point to the same raw data as the original.
  my $trf_dataset = $dataset->copy ( { name => undef } );
  $trf_dataset->data ( $dataset->data() );

  EP::Core::Data::Parser::exp::user->_write_dataset_data ( $trf_dataset, $trf_dataset->data );
  my $meta = $trf_dataset->data->process_metadata();

  # HACK. This is to force any future read of raw data on the original (pre-selection) dataset to re-read it from the
  # disk, rather than read the modified raw data from trf_dataset.
  delete $dataset->{data};

  $trf_dataset->data->raw_data->badflag(1);
  return ($meta, $trf_dataset);
}

# sub: transpose_data
#   Changes rows and columns of the expression matrix.
sub transpose_data {
  my $self = shift;
  my ($dataset) = @_;

  my $matrix = $dataset->data->raw_data;

  my $tmp_dataset_row   = $dataset->data->{row};
  $dataset->data->{row} = $dataset->data->{column};
  $dataset->data->{column} = $tmp_dataset_row;
  $dataset->data->{matrix} = transpose ( $matrix );
}

# sub: impute_data_knn
#   Performs K-nearest neighbours data imputation via R.
sub impute_data_knn {
  my $self = shift;
  my ($dataset, $nbrs) = @_;

  my $matrix = $dataset->data->raw_data;

  my ($folder) = $EP::Config::h->{EPNG_USRDATA_PATH} . "/" . $dataset->folder->location;
  my ($filename) = $dataset->filename;

  my $dataset_filepath = $folder . "/" . $filename;
  my $tmp = new File::Temp;

  EP::Rclient->R(<<__RCODE__);
app_root = "$EP::Config::h->{ROOT}";
source("$EP::Config::h->{ROOT}/cgi-lib/R/ep.io.R");
library(pamr);
M = ep.readBin("$dataset_filepath.bin");
O = pamr.knnimpute( list( x = M ), k = $nbrs );
ep.writeBin( O\$x, "$tmp" );
__RCODE__

  my ( $r, $c, $data ) = readflex ( "$tmp" );
  warn "KNN IMPUTED: ", $data->info if DEBUG;
  $dataset->data->{matrix} = $data;
}

# sub: impute_data_replace_with_rowavgs
#   Replaces NAs with row averages
sub impute_data_replace_with_rowavgs {
  my $self = shift;
  my ($dataset) = @_;

  my $matrix = $dataset->data->raw_data;
  $dataset->data->{matrix}->where($matrix->isbad) .=
    (transpose($matrix->average) x ones($matrix->dim(0)))->where($matrix->isbad);
}

# sub: impute_data_replace_with_zeroes
#   Replaces NAs with zeroes.
sub impute_data_replace_with_zeroes {
  my $self = shift;
  my ($dataset) = @_;

  # TODO: CHECK THIS (setnantoval?)
  my $matrix = $dataset->data->raw_data;
  $dataset->data->{matrix}->inplace->setbadtoval(0);
}

# sub: avg_row_id
#   Averages values of duplicate row identifiers.
sub avg_row_id {
  my $self = shift;
  my ($dataset) = @_;

  my $matrix = $dataset->data->raw_data;
  my ( %uniq_rows, $data, @matrix, $rows );

  my @row_annot_names = $dataset->data->{row}{0}{full_annot}->Keys;

  #get the data for each unique identifier
  my $row_count = 0;
  foreach my $row ( sort { $a <=> $b } keys %{ $dataset->data->{row} } ) {
    my $header = $dataset->data->{row}{$row}{full_annot}->Values ( $row_annot_names[0] );
    push @{$uniq_rows{$header}}, $row;
  }

  my $d = transpose ( $matrix );
  my @new_data;
  $row_count = 0;
  foreach my $row_header ( sort keys %uniq_rows ) {
    my $row_idx = pdl @{$uniq_rows{$row_header}};
    push @new_data, $d($row_idx,:)->average;

    $rows->{$row_count}{full_annot} = $dataset->data->{row}{$row_idx->at(0)}{full_annot};
    $rows->{$row_count}{header}     = $dataset->data->{row}{$row_idx->at(0)}{header};
    $row_count++;
  }

  #remove old data and put in new
  $dataset->data->{row} = $rows;
  $dataset->data->{matrix} = cat @new_data;
}

# sub: ratio_to_logratio
#   Converts matrix of ratios to matrix of log N ratios.
sub ratio_to_logratio {
  my $self = shift;
  my ($dataset, $log) = @_;

  my $matrix = $dataset->data->raw_data;

  $dataset->data->{matrix}->inplace->log();
  $dataset->data->{matrix}->inplace->divide ( log ( $log ), 0 ) unless $log eq "e";
  $dataset->data->{matrix}->inplace->setnantobad();

  my @col_annot_names = $dataset->data->{column}{0}{full_annot}->Keys;

  foreach my $column ( sort { $a <=> $b } keys %{ $dataset->data->{column} } ) {
    $dataset->data->{column}{$column}{full_annot}->Replace ( 0 , "Log $log " . "[" .  $dataset->data->{column}{$column}{full_annot}->Values ( 0 ) . "]" );
  }

  return $dataset;
}

# sub: intensity_to_ratio
#   Converts matrix of intensity data to matrix of ratios (log or normal).
sub intensity_to_ratio {
  my $self = shift;
  my ($dataset, $nr_columns, $ch1, $ch2, $log, $dye_swaps) = @_;

  my $matrix = $dataset->data->raw_data;

  my (@matrix, $column);

  die "Insufficient parameters for transformation" if !$nr_columns || !$ch1 || !$ch2;

  my $nr_exp = scalar( keys %{ $dataset->data->{column} } ) / $nr_columns;
  die "Number of experiments do not match with number of columns in data matrix" if $nr_exp != int ($nr_exp);

  #get the dye swap exp
  my @ds;
  foreach my $dye_swap ( split /,/, $dye_swaps ) {
    $dye_swap =~ s/^\s+|\s+$//g;
    next if $dye_swap == 0 && $dye_swap ne "0"; # EP::Common::General::data_type ( $dye_swap ) != NUMERIC;
    push @ds, ($dye_swap - 1);
  }

  my $exp_idx = sequence ( $nr_exp ) * $nr_columns;
  my $ch1_idx = $exp_idx + $ch1 - 1;
  my $ch2_idx = $exp_idx + $ch2 - 1;
  my $swp_idx = pdl (@ds);
  
  warn "ch1: $ch1_idx\nch2: $ch2_idx";

  my ($power) = $log =~ /none_(.*)/;

  $matrix = &pow ( $power, $matrix->inplace ) if defined $power and $power ne "e";
  $matrix->inplace->exp if defined $power and $power eq "e";

  warn "Swapping $swp_idx";
  my $r = $matrix($ch1_idx,:) / $matrix($ch2_idx,:);
  $r($swp_idx,:) .= 1/$r($swp_idx,:);

  if ( $log ne 'none' ) {
    $log = $power if $log =~ /none/;
    $r->inplace->log();
    $r->inplace->divide ( log ( $log ), 0 ) unless $log eq "e";
    $r->inplace->setnantobad;
  }

  $r->inplace->setnantobad();

  #set the new column headers
  foreach ( 0 .. $exp_idx->dim(0) - 1 ) {
    my $exp     = $exp_idx->at($_);
    my $ch1_col = $ch1_idx->at($_);
    my $ch2_col = $ch2_idx->at($_);

    $column->{$exp}{full_annot} = $dataset->data->{column}{$ch1_col}{full_annot}; # assume annotation is the same for both ?
    $column->{$exp}{full_annot}->Replace ( 0, 'Log ' . $log . ' [' . $dataset->data->{column}{$ch1_col}{full_annot}->Values ( 0 ) . " / " . $dataset->data->{column}{$ch2_col}{full_annot}->Values ( 0 ) . ']' );
  }

  #remove old data and put in new
  $dataset->data->{column} = $column;
  $dataset->data->{matrix} = $r;
}

# sub: affy_rel_ratio
#   Converts Affy absolute values to relative ratios, dividing each gene's value by its overall mean across all
#   conditions. Optionally can also take log of these values.
sub affy_rel_ratio {
  my $self = shift;
  my ($dataset, $log, $rel_to, $rel_to_column_name) = @_;

  my @rel_to_column_names = split /\r\n|\r|\n/, $rel_to_column_name;

  my $matrix = $dataset->data->raw_data;
  my @rel_to_column;

  for my $col ( 1 .. $matrix->dim(0) ) {
    push @rel_to_column, $col if grep { $dataset->data->{column}{$col-1}->{full_annot}->Values ( 0 ) eq $_ } @rel_to_column_names;
  }

warn "Looking for @rel_to_column";

  my ($power) = $log =~ /none_(.*)/;

  if ( $rel_to eq 'gene_mean' ) {
    warn "First we will unlog this if necessary... power = $power";
    if ( defined $power ) {
      $matrix->inplace->mult ( log ( $power ), 0 ) if $power ne "e";
      $matrix->inplace->exp;
    }
    $matrix->inplace->divide ( transpose ( $matrix->average ), 0 );
  } elsif ( $rel_to eq 'column' ) {
    die "Need to specify valid column number (between 1 and " . $matrix->dim(0) . ")!" if $#rel_to_column == -1 or any pdl (@rel_to_column) > $matrix->dim(0);

    my $s = -1 + pdl(@rel_to_column);
    my $m = (defined $power) ? pow ( $power, $matrix($s,) ) : $matrix($s,);
    die "Transformation failed. Maybe these data are not already log-transformed?" if any not isfinite $m;
    my $d = $m->sumover->transpose->divide ( scalar @rel_to_column, 0 );
    
    if ( not defined $power ) {
      $matrix->inplace->divide ( $d, 0 );
    } else {
      $d->inplace->log();
      $d->inplace->divide ( log ( $power ), 0 );
      $matrix->inplace->minus ( $d, 0 );
    }

    if ( scalar @rel_to_column == 1 ) {
      for my $col ( 1 .. $matrix->dim(0) ) {
        $dataset->data->{column}{$col-1}{full_annot}->Replace ( 0, $dataset->data->{column}{$col-1}{full_annot}->Values ( 0 ) . " rel to " . $dataset->data->{column}{$rel_to_column[0]-1}{full_annot}->Values ( 0 ) ) if $col != $rel_to_column[0];
      }
    }

    delete $dataset->data->{column}{$_-1} foreach @rel_to_column;

    my $q = which(not(sequence($matrix->dim(0))->in(-1 + pdl @rel_to_column)));
    $dataset->data->{matrix} = $matrix($q,);

    return if defined $power;
  }

  if ( $log ne "none" ) {
    $log = $power if $log =~ /none/;

    warn "Take log ($log) of the data now";
    $matrix->inplace->log();
    $matrix->inplace->divide ( log ( $log ), 0 );
    $matrix->inplace->setnantobad;
  }
}

# sub: mean_center
#   Converts Affy absolute values to relative ratios, dividing each gene's value by its overall mean across all
#   conditions. Optionally can also take log of these values.
sub mean_center {
  my $self = shift;
  my ($dataset, $mean_center) = @_;

  my $matrix = $dataset->data->raw_data;

  $matrix .= $matrix - $matrix->average->transpose if $mean_center =~ /genes/;
  $matrix .= $matrix - $matrix->transpose->average if $mean_center =~ /arrays/;
}

1;
